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ABSTRACT 

Supernova (SN)-driven pregalactic outflows may be an efficient mechanism for distributing the product 
of stellar nucleosynthesis over large cosmological volumes prior to the reionization epoch. Here we 
present results from three-dimensional numerical simulations of the dynamics of SN-driven bubbles as 
they propagate through and escape the grasp of subgalactic halos with masses M — 10 8 hr 1 M Q at 
redshift z = 9. Halos in this mass range are characterized by very short dynamical timescales (and even 
shorter gas cooling times) and may therefore form stars in a rapid but intense burst before SN 'feedback' 
quenches further star formation. The hydrodynamic simulations use a nested grid method to follow the 
evolution of explosive multi-SN events operating on the characteristic timescale of a few xl0 7 yr, the 
lifetime of massive stars. The results confirm that, if the star formation efficiency of subgalactic halos is 
10%, a significant fraction of the halo gas will be lifted out of the potential well ('blow-away'), shock 
the intergalactic medium, and pollute it with metal enriched material, a scenario recently advocated by 
Madau, Ferrara, & Rees (2001). The volume filling factor of the ejecta is of order unity. Depending on 
the stellar distribution, we find that less than 30% of the available SN energy gets converted into kinetic 
energy of the blown away material, the remainder being radiated away. It appears that mechanical 
feedback is less efficient than expected from simple energetic arguments, as off-nuclear SN explosions 
drive inward-propagating shocks that tend to collect and pile up cold gas in the central regions of the 
host halo. Low-mass galaxies at early epochs may survive multiple SN events and continue forming 
stars. 

Subject headings: cosmology: theory - galaxies: formation - hydrodynamics - intergalactic medium - 
supernovae: general 
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1. INTRODUCTION 

In currently popular hierarchical clustering scenarios for 
the formation of cosmic structures - all variants of the 
cold dark matter (CDM) cosmogony - the assembly of 
galaxies is a bottom-up process in which large systems 
result from the merging of smaller subunits. In these the- 
ories dark matter halos with masses M = 10 s h~ 1 M.Q 
collapse at z » 10 from 2-cr density fluctuations: the 
gas infalling along with the dark matter perturbation is 
shock-heated to the virial temperature T vu = 2 x 10 4 K, 
condenses rapidly due to atomic line cooling, and becomes 
self-gravitating. Massive stars then form with some ini- 
tial mass function (IMF), synthesize heavy-elements, and 
explode as supernovae (SNe) after a few x 10 7 yr, enrich- 
ing the surrounding medium. These subgalactic stellar 
systems, possibly aided by a population of accreting black 
holes in their nuclei and/or by an earlier generation of stars 
in even smaller halos ('minihalos' with virial temperatures 
of only a few hundred kelvins, where collisional excitation 
of molecular hydrogen is the main coolant), are believed 
to have generated the ultraviolet radiation and mechanical 
energy that reheated and reionized the universe (see Loeb 
& Barkana 2001 for a recent review). 

It is a simple expectation of the above scenario that the 
energy deposition by SNe in the shallow potential wells of 
subgalactic systems may have two main effects, depending 



on the efficiency with which halo gas can cool and fragment 
into clouds and then into massive stars: (i) the disrup- 
tion of the newly formed object, i.e. the most violent ver- 
sion of the so-called 'stellar feedback' (Larson 1974; Dekel 
& Silk 1986; Mori et al. 1997; Mori, Yoshii, & Nomoto 
1999; MacLow & Ferrara 1999; D'Ercole & Brighenti 
1999; Murakami & Babul 1999; Ferrara & Tolstoy 2000; 
Ciardi et al. 2000); and (ii) the blow-away of metal- 
enriched baryons from the host (dwarf) galaxy, causing 
the pollution of the IGM at early times (Tegmark, Silk, & 
Evrard 1993; Voit 1996; Nath & Trentham 1997; Gnedin 
& Ostriker 1997; Ferrara, Pettini, & Shchekinov 2000; 
Madau, Ferrara, & Rees 2001, hereafter Paper I). The 
well-established existence of heavy elements like carbon, 
nitrogen, and silicon in the Lya forest clouds at z = 3 — 3.5 
may be the best evidence for such an early episode of pre- 
galactic star formation and outflows (Songaila 1997; El- 
lison et al. 2000). Stellar feedback and galactic outflows 
at high-redshifts may also temporarily halt or delay the 
formation of dwarf galaxies (Scannapieco & Broadhurst 
2001; Scannapieco, Ferrara, & Broadhurst 2000), affect 
the thermal state of the IGM (Madau 2000; Theuns, Mo, 
& Schaye 2001; Cen & Bryan 2001), and play a crucial role 
during the reionization epoch (Ciardi et al. 2000; Bruscoli 
et al. 2000). 

This is the second paper in a series aimed at a detailed 
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study of stellar feedback and SN-driven pregalactic out- 
flows at early epochs. In Paper I it was shown that the ob- 
served narrow Doppler widths and inferred large filling fac- 
tor of chemically enriched, low-density Lya forest clouds 
may point to a more uniform (i.e. 'early') rather than in 
situ (i.e. 'late') cosmic metal pollution. Whilst outflows 
of metal-rich gas are directly observed in Lyman-break 
galaxies at z w 3 (Pettini et al. 2000), most of this gas may 
not leave the immediate surroundings of these deep grav- 
itational potential wells. We argued that, if intergalactic 
metals were actually dispersed over large cosmological vol- 
umes from massive galaxy halos at late times, such a de- 
layed epoch of galactic "super-winds" would have severely 
perturbed the IGM (since the kinetic energy of the ejecta is 
absorbed by intergalactic gas) , raising it to a much higher 
adiabat than expected from photoionization, and produc- 
ing large spatial variations of the baryons relative to the 
dark matter. This would make the success of cosmolog- 
ical hydrodynamical simulations in matching the overall 
observed properties of Lya absorption systems largely co- 
incidental. Alternatively, the observations could be best 
explained if the ejection of heavy elements at velocities 
exceeding the small escape speed of numerous subgalactic 
systems were to take place at very high redshifts. This 
is because hot enriched material cools more efficiently at 
these early epochs, and the expansion of SN-driven bub- 
bles into a dense IGM (pre-photoionized by the same mas- 
sive stars that later explode as SNe) would be halted by 
the external pressure. In this scenario, any residual pecu- 
liar velocity would have been rcdshiftcd away by z = 3, 
the Lya forest would be hydrodynamically 'cold', and the 
intergalactic baryons would have relaxed again under the 
influence of dark matter gravity. 

In Paper I the most efficient pollutant of the IGM on 
large scales were identified as subgalactic systems with 
masses ~ 10 8 M Q at z <~ 10, when large number of 
them grow non-linear and collapse. We showed how stellar 
feedback in these pregalactic halos (with virial tempera- 
tures T v ; r = 2 x 10 4 K, i.e. at the peak of the cooling curve 
for primordial gas) may be an inherently different pro- 
cess than that often invoked to prevent the 'overcooling' of 
baryons in the central regions of larger systems (White & 
Frenk 1991; Cole et al. 1994; Navarro & Steinmetz 1997, 
2000). This is because such fragments are characterized 
by very short dynamical timescales (and even shorter gas 
cooling times): throughout the halo the accreted baryons 
condense immediately due to atomic hydrogen cooling, 
and the supply of cold gas for star formation is only lim- 
ited by the infall rate. They are therefore expected to go 
through a rapid but intense star forming phase, since one 
cannot assume instantaneous feedback from SNe in systems 
with such a short dynamical timescalc. Lacking the ability 
to store efficiently the corresponding energy input in tur- 
bulent and magnetic forms or by relaxing to a multiphase 
medium, these system have very limited resources to sus- 
tain a self-regulated star formation cycle (cf. Efstathiou 
2000). 

In Paper I we used the thin shell approximation to 
model the evolution of SN-driven bubbles as they propa- 
gate through and blow out from subgalactic halos. We 
assumed spatial coherency among SN events and tem- 

5 Unless otherwise stated, we will assume throughout this paper an 
constant of Hq = 100 h km s _1 Mpc -1 . The numerical simulations are 



poral coherency among their progenitors, i.e. all explo- 
sions were assumed to take place at the center of the 
halo, and the progenitor stellar population to form co- 
evally, on a timescale short compared to the lifetime of 
massive stars. In this paper we present results from three- 
dimensional hydrodynamic simulations of the dynamics of 
SN-driven bubbles in subgalactic halos using a hierarchical 
nested grid technique. Besides accounting for the mass- 
dependent main sequence evolutionary timescale of mas- 
sive stars (as in Paper I), which causes a temporal spread 
of SN explosions, we also explore the role of different spa- 
tial distributions of the explosion sites inside the galaxy, 
thus relaxing the assumption that all SNe occur at the 
center. The crucial questions we try to give a quantitative 
answer in this paper are: How efficient mechanical feed- 
back really is? Can small halos survive multiple SN events 
and continue forming stars or are they blown away and 
end their life as 'naked stellar clusters' ? How far from the 
production sites can metals be ejected into the IGM? 

2. MODEL ASSUMPTIONS 

We will model the structural properties of subgalactic 
systems following Paper I. We shall assume that virial- 
ized dark matter halos, formed through hierarchical clus- 
tering, have a universal (spherically averaged) Navarro- 
Frenk-White (1997, hereafter NFW) density profile, 
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where x = r/r vir , r vir is the virial radius of the sys- 
tem, i.e. the radius of the sphere encompassing a mean 
overdensity of 200, c is the halo concentration parame- 
ter, S c — (200/3)c 3 /F(c) is a characteristic overdensity, 
Pcrit = 3iJ 2 / (8ttG) the critical density of the universe 5 at 
rcdshift z, and 
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The total within the virial radius 
(47r/3)200p cr itr^ ir , and the 'typical' concentration param- 
eter of halos at z = 10 with M = 10 8 /i _1 M Q is c = 5 
(Paper I). Note that high-resolution N-body simulations 
by Bullock et al. (2001) indicates that high-redshift ha- 
los are actually less concentrated that expected from the 
NFW prediction. In this case we may be overestimating 
the escape speed from subgalactic systems. 

For our prototypical pollutor (below the mass scale given 
above the main coolant is molecular hydrogen and halos 
are strongly affected by radiative feedback, while above 
this mass scale cooling is less rapid, halos are more rare, 
and metals are more efficiently retained by the deeper po- 
tential well, see detailed discussion in Paper I), the charac- 
teristic virial radius, virial temperature, circular velocity 
at the virial radius, and escape velocity at the center are 



= 0.75 kpcM 8 1/3 /i 



2 x 10 4 KM. 
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Einstcin-de Sitter cosmology with Qm = 1, S7a = 0, and a Hubble 
run with h = 0.5. 
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respectively, where M§ is the halo mass in units ol 
10 8 h^ 1 M©. Here we have set the mean molecular weight 
/U to 0.59, appropriate for a fully ionized primordial gas. 
At the virial radius v e (r v i r ) = 0.62 v e (0). If baryons virial- 
ize in the dark matter halo to an isothermal distribution, 
they will be shock-heated to the virial temperature and 
settle down to a density profile 
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where m p is the proton mass. The central gas density po is 
determined by the condition that the total baryonic mass 
fraction within the virial radius is equal to initially. 
Adopting n b h 2 = 0.019 (Buries & Tytler 1998), one gets 
po = 840/i~ 2 p crit = 1.6 x 10~ 23 gr cm" 3 at z = 9. The 
halo gas density at r vir is p(r v i r ) = 0.00144p = 2.3 x 10~ 26 
gr cm -3 . As this is about 60 times higher at (1 + z) = 10 
than the IGM density, to avoid unphysical effects due to 
this jump we have allowed the two distributions to merge 
trough a hyperbolic tangential transition of width 0.05 r v ; r . 

Numerical N-body/hydrodynamics simulations of struc- 
ture formation in CDM cosmologies have provided in the 
last few years a definite picture for the topology of the 
IGM, one of an interconnected network of sheets and fila- 
ments with virialized systems located at their points of 
intersection. Intergalactic gas will be infalling onto a 
galaxy along the filaments of this "cosmic web", at ap- 
proximately the halo escape velocity. Our simulations 
do not include the effect of this non-spherical infall on 
SN-driven outflows. Note that, for the halos/outfiows un- 
der consideration, the escape speed at the virial radius is 
v e (r v ir) = 48 kms -1 , while the outflow velocity is close to 
150 kms -1 for a star formation efficiency of 10% (see Fig. 
7 of Paper I). In this case the pressure due to infalling 
material is 9 times smaller than the ram pressure of the 
expanding shell, and can be neglected in first approxima- 
tion. As the blast propagates into the IGM, however, it 
may preferentially enrich the voids in between the denser 
filaments. SN-driven outflows in a realistic cosmological 
density field will be the subject of a subsequent paper in 
this series. 

2.1. Supernova progenitors 

As shown in Paper I, for halo masses in the range 
10 8 h' 1 < M < 10 10 h' 1 M Q at (1 + z) = 10, the cooling 
time is shorter than the dynamical timescale everywhere 
in the halo: infalling gas never comes to hydrostatic equi- 
librium, but collapses to the center at the free-fall rate. 
In our prototypical halo, the free-fall time is 10 7 yr at 
r = 0.1r v ; r , increasing to 8 x 10 7 yr at r v ; r (Paper I). 
The cooling time for gas at T w -„ is always shorter than the 
free-fall time by more than two orders of magnitude. Me- 
chanical energy will be injected by SNe only after a few 



times 10 yr: at this stage SN-driven bubbles will prop- 
agate into the halo quenching further star formation, and 
the conversion of cold gas into stars will be limited by the 
increasing fractional volume occupied by supernova rem- 
nants. 

Before SN feedback starts operating, some fraction /* 
of the halo gas may be able to cool, fragment, and form 
stars. In general, cooling is efficient when the cooling time 
£ C ooi = 1.5nfeT/(nijA) is much shorter than the local gas- 
dynamical timescale idyn = ('InGp)^ 1 ^ 2 , where A is the 
radiative cooling funtion. This condition implies that the 
energy deposited by gravitational contraction cannot bal- 
ance radiative losses; as a consequence, the temperature 
decreases with increasing density, the cloud cools and then 
fragments. At any given time, fragments form on a scale 
that is small enough to ensure pressure equilibrium at the 
corresponding temperature, i.e. the Jeans length scale. 

Small pregalactic systems at early times, like the ones 
under consideration, may then be expected to consume 
a fair fraction of their cold gas in a single burst of star 
formation. Note that the highest efficiency estimated in 
nearby star-forming regions is about 30% for the Ophiuchi 
dark cloud (Wilking & Lada 1983; Lada & Wilking 1984), 
and that the UV radiation from massive stars may further 
inhibit cooling and delay the collapse of mildly overdense 
regions (Nishi & Tashiro 2000). In the following we will 
adopt an illustrative value of /* = 10% for the star for- 
mation efficiency. This is consistent with the stellar-to- 
baryonic mass ratio observed today, and implies that only 
a few percent of the present-day stellar mass would form 
at these early epochs (see discussion in §6). The ensu- 
ing SNe (assuming a Salpeter initial mass function) would 
pollute the IGM to a mean metallicity of order 0.3% solar, 
comparable to the levels observed in the Lya forest at red- 
shift 3. Given our poor understanding of star formation, 
however, /* should really be considered at this stage as a 
free parameter of the theory. 

The halo under study has a baryonic mass M& = Q&M = 
0.019 h- 2 x 10 8 /i -1 M Q = 1.9 x 10 6 h- 3 M o . Its total stel- 
lar mass is then M* = f±M b = 1.5 x 10 6 M Q (h = 0.5). 
An amount M* of baryons is subtracted from the total gas 
mass; the rest of the gas is assumed to keep the same initial 
distribution. We assume that the initial stellar mass func- 
tion (IMF) can be approximated by a Salpeter power-law 
with lower and upper mass cutoff's equal to mi — 0.1 M© 
and m u = 120 M©, respectively. Note that considerable 
uncertainties still remain on the characteristic mass of 
the first luminous objects. Numerical simulations of the 
fragmentation of primordial clouds with masses ^ a few 
x 10 5 Mq in hierarchical cosmogonies have suggested that 
the IMF of the very first, zero-metallicity stars forming at 
z £ 20 may be extremely top-heavy (Bromm, Coppi, & 
Larson 1999, 2001a; Abel, Bryan, & Norman 2000), per- 
haps giving origins to a population of pregalactic massive 
black holes (Madau & Rees 2001). As this feature appears 
linked to primordial H 2 chemistry and cooling, at the later 
epochs of interest here and in more massive halos where 
gas condense due to atomic line cooling, it is plausible that 
a 'second' generation of stars may be able to form with an 
IMF that is less biased towards very high stellar masses 
(Bromm et al. 2001b). 

In our simulation, the mass of each star is determined 
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by randomly sampling the IMF until the total stellar mass 
is equal to M*; SN progenitors are then identified as those 
stars more massive than 8 Mq . We find that a total num- 
ber of N t = 11170 SNe are produced in the system, out of 
4.24 x 10 6 stars (average stellar mass = 0.35 M Q ) yield- 
ing a SN every 134 Mq of stars formed. Each massive 
star (with mass m in units of Mq) explodes after a main 
sequence lifetime of 

log^ms (yr) = 10.025 - 3.559 log m + 0.898(logm) 2 . (8) 

This is a fit to a compilation of the available data in the 
literature (Schaller ct al. 1992; Vacca, Garmany & Shull 
1996; Schaerer & de Koter 1997; Palla, private communica- 
tion) . In the local universe SNe are known not to occur at 
random location but rather to cluster into OB associations 
(Heiles 1990). In nearby galaxies the luminosity function 
of OB associations is well approximated by a power-law 
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with (3 « 2 (McKee & Williams 1997; Oey & Clarke 1998). 
Here Nob is the number of associations containing TV OB 
stars; the probability for a cluster of OB stars to host 
N SNe is then cx N~ 2 . We apply a grouping procedure 
according to the above distribution: this yields 90 OB as- 
sociations, each containing a variable number of massive 
stars ranging from a few tens to up to 2000. For numerical 
reasons to be discussed below, we set a lower limit to the 
number of stars in an association equal to 20. 

As a final step we have to spatially locate the 90 asso- 
ciations in the halo. Our modelling of star formation is, 
unfortunately, limited by finite numerical resolution and 
the neglect of self-gravity. We have therefore decided to 
use a simple scheme that incorporates some gross observed 
features of star formation and makes use of actual gas 
properties such as the local mass density. We distribute 
SNe according to the gas density power-law introduced 
by Schmidt (1959), i.e. the number density of explosions 
is proportional to p a . In starburst and spiral galaxies, the 
disk-averaged star formation rates and gas densities are 
well represented by a Schmidt law with a — 1.4 ± 0.15 
(Kennicutt 1998). Numerous theoretical scenarios that 
produce a Schmidt law with a = 1 — 2 can be found in the 
literature (e.g. Silk 1997); since the cooling time is much 
shorter than the gas dynamical timescale in dense star- 
forming regions, the star formation rate is often taken to 
be proportional to p/idyn oc p 3 / 2 . Similar phenomenologi- 
cal prescriptions for star formation are used in numerical 
simulations of hierarchical galaxy formation (Katz 1992). 
In this paper we run two simulations, corresponding to a 
relatively extended (a = 1) or more centrally concentrated 
(a = 2) star formation volume (in the limit a — > oo all SNe 
explode at the center). We will refer to them in the fol- 
lowing as Case 1 (a = 1) and Case 2 (a = 2), respectively. 

It is important to notice that our assumption of an ex- 
tended spherical distribution of SNe may be a poor one if 
the fragmentation of cooling gas freely-falling towards the 
halo center is an inefficient process (as argued by Kash- 
linsky & Rees 1983), and star formation only occurs af- 
ter cold halo gas actually settles down in a centrifugally- 
supported galactic disk. If an exponential disk with scale 
length rd forms, and the specific angular momentum of the 



disk material is the same as that of the dark halo (treated 
for simplicity as a singular isothermal sphere), then an- 
gular momentum conservation fixes the collapse factor to 
fvirAd = \/2/A, where A rss 0.05 is the spin parameter. 
Our fiducial halo/disk system at z — 9 would then be char- 
acterized by a scale length of only r d = 27/i _1 pc (Paper 
I). In this scenario, and with the resolution of our numeri- 
cal simulations (see below) , all SNe explosions would take 
place within a few central mesh cells. Also, a thin dense 
disk would hardly couple to the outflow, and matter would 
be ejected perpendicularly to the disk (MacLow & Ferrara 
1999). It is unclear, however, whether thin self-gravitating 
disks would actually form and/or survive in subgalactic 
fragments. If a disk of mass M d follows an isothermal ver- 
tical profile with a thermal speed c s = 10 kms , typical 
of gas which is continuosly photohcatcd by stars embed- 
ded within the d isk itself, then its scale height at radius r d 
is h/r d = Vl6eX(M/M d )(c s /V c ) 2 , and the disk would be 
thick rather than thin. Recent numerical simulations of 
the formation and fragmentation of primordial molecular 
clouds in CDM cosmologies, with realistic initial condi- 
tions, do not show the presence of a disk, and form the 
first dense fragments close to the center of the halo (Abel 
et al. 2000). A rotationally-supported disk that fragments 
efficiently appears to only form as a consequence of ideal- 
ized initial conditions, i.e. a top-hat isolated sphere that 
initially rotate as a solid body (Bromm et al. 1999). 

3. PRELIMINARY CONSIDERATIONS 

Before discussing the detailed results of the numerical 
simulations, it is important to provide some physical in- 
sights into the problem. 

3.1. Effects of explosion location 

In order for the gas to be ejected from the galaxy, a nec- 
essary condition implies that the initial energy provided by 
the supernovae has to be larger than the gravitational one. 
The latter is given by 



(10) 



where </>(r) = —GM(r)/r is the gravitational potential. 
With the gas density distribution given by equation (R), 
integration yields 
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This energy must be compared with the total mechan- 
ical energy injected by SN explosions, = N t E = 
1.1 x 10 55 erg. The naive conclusion would be that this en- 
ergy deposition, being roughly hundred times higher than 
the gas binding energy, should produce a complete blow- 
away of the halo. There are two main reasons, however, 
for why this is not actually the case: i) the conversion 
efficiency of E^n m to kinetic energy of the interstellar 
medium (ISM) is well below unity since radiative energy 
losses, particularly at the halo center where the baryon 
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density is initially n(0) w 10cm~ 3 , carry away a signifi- 
cant fraction of this energy; and ii) as the explosion sites 
are scattered within the host galaxy, they can have very 
different effects. For example, explosions taking place in 
the outer halo will tend to expand outwards along the 
rapidly decreasing density gradient: the released energy 
will be eventually used to energize the 1GM. Off-nuclear 
SN explosions will drive inward-propagating shocks that 
pile up and compress the gas in the central regions, fur- 
ther increasing radiative losses. Also, energy deposition by 
an association at the center may never accelerate material 
beyond the virial radius. 

We can better quantify the last point by adopting the 
following zeroth order approximation. We treat each asso- 
ciation as a point explosion (in practice, due to the spread 
of the various SN events, the ensuing superbubble is better 
described by a wind solution, see Paper 1). The evolution of 
a point explosion in a stratified medium can be treated by 
the thin shell approximation. This solution accurately ap- 
proximates the exact numerical results (see, e.g., MacLow 
et al. 1989) and can be obtained from dimensional anal- 
ysis. The shell velocity is v s (r) ~ [P/p(r)] 1 ' 2 , where the 
pressure P is roughly equal to E/X 3 : E = NE is the 
energy injected by N SNe in the association. The explo- 
sion is allowed to take place at different radial positions 
x o = r o Avir ; the bubble size along the radial coordinate is 
then X = r V i T (x — xq). We can then calculate the value of 
«s(?"vir)) which can then be compared to the escape veloc- 
ity u e (r v ; r ) to get an idea of the final fate (escape versus 
recapture) of the ejected gas. The results are shown in 
Figures @ and | as u s (Tvir) isocontours in the plane iden- 
tified by xq and N. Also shown are the locations in that 
plane of the associations used in the numerical simulations 
for the two cases considered. It can be seen that, if the 
explosion occurs at the halo center (xq = 0), a minimum 
of about 50 SNe is required in order to achieve shell veloci- 
ties in excess of v e (r v i T ). This constraint is much weaker at 
(say) xo — 0.7, where even a single SN event can produce 
the same effect. Thus, every explosion occurring beyond 
this location will easily lead to mass loss from the galaxy. 
Note that in Case 1 a larger fraction (about 90%) of SN 
events will give origin to significant mass loss, as they take 
place preferentially in the galactic outskirts. By contrast, 
in Case 2 about 1 /3 of the OB associations generate only 
sub-escape speeds, i.e. v a (r v i T ) < v e (r v - n -). 

3.2. Energy- and mass-carrying winds 

Explosions located in the halo outer regions produce 
higher velocities of the outflowing material. The swept- 
up mass, however, decreases as the halo density is lower 
in the outskirts. In this case SNe basically vent energy 
(an energy-carrying wind) into the IGM, but little galac- 
tic mass. An accurate determination of the kinetic energy 
carried away, Ek oc M s v 2 , can only be done numerically as 
it depends on the detailed mass loading of the propagating 
shells. To get some insight, we estimate the kinetic energy 
flux through r v i r as 

£ k (x ,N) = v 2 s (r vi rMx Q ) = / dxp, (12) 

A J p J Xo 

where S(xo) is the halo gas (mass) column density along 
the radial direction from xq to r v j r . We then plot in Figure 



|3| the column £ normalized to the total column through 
the halo center, E t , as a function of the square of the ratio 
Vs ( r vir ) / v e (»Vir ) ■ If the kinetic energy flux were the same at 
every location one would find Sw 2 (r V i r ) = const = £k(N). 
As seen from Figure however, the trend is rather dif- 
ferent. For N — 10 the curve drops very rapidly below 
the one for £^(10): this is because the increase in the shell 
velocity cannot compensate for the drop in mass as the 
explosion is displaced from the origin. It is only in the 
outskirts (xq > 0.75) that the transition to a predomi- 
nantly energy-carrying wind occurs. In the N — 100 case 
such inversion in not seen, essentially because the mass 
decrease always dominate over the velocity increase. 

In conclusion, if SNe occur in small N ss 20 — 30 asso- 
ciations (or are isolated), either there is no mass-loss (if 
they are located at the center) or they mostly inject en- 
ergy rather than mass into the IGM (if they are located in 
the outer halo). It is only when they are grouped together 
in larger associations that they will be able to eject both 
mass and energy into the intergalactic space. 

3.3. Porosity and amplification effect 

Based on the results shown in Figures |l| and ||, Case 
2 should lead to a smaller mass loss from our subgalac- 
tic halo. This, however, may not be necessarily the case, 
as "amplification" effects can be important in the central 
regions. Due the high stellar density at x < 0.2 (say), 
overlapping of the hot bubble interiors does occur if cool- 
ing is not too strong. In this case, the energy deposited at 
the various explosion sites can add and act coherently on 
the same shell, leading to an amplification of the pressure 
force and to higher shell velocities. In other words, the 
porosity of the hot gas, Q, (its volume filling factor being 
equal to 1 — e~®) is a key parameter, as already pointed 
out by Silk (1997) and Efstahtiou (2000) in the context 
of feedback-regulated star formation in galaxy disks. The 
effect of the interaction and merging of individual bub- 
bles will be discussed in the framework of the numerical 
simulations discussed below. 

4. NUMERICAL METHOD 

The evolution of the gas is described by the three 
dimensional hydrodynamic equations for a perfect fluid 
in Cartesian geometry. The continuity equation (includ- 
ing a term due to mass ejection from SNe, p s ), the mo- 
mentum equation with the gravitational acceleration g = 
(9x, 9y: gz) T > and the thermal energy equation associated 
with the rates of cooling A and SNe heating T s , can be 
written in compact form as 



with 





dV dF dG OH 

dt dx dy dz 


(13) 


u = 


(p, pv x , pv y , pv z , pe) T ', 


(14) 


F = 


(pv x , pv 2 x + p, pv x v y ,pv x v z ,pv x h e ) T , 


(15) 


G = 


(pVy,pV X Vy,pvl pVyV Z ,pVyh e ) T , 


(16) 


H = 


[pV z , pv x v z , pv y v z , pv\ + p, pv z h e ) T , 


(17) 


S = 


(p s , P9x, pg y , pg z , pv ■ g + r s - A) T . 


(18) 
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Here U is a state vector of conserved quantities, F, G, 
and H are the corresponding flux vectors, and S is the 
source-term vector that includes sources and sinks of con- 
served quantities, such as heating and cooling terms and 
the gravitational acceleration. Also, p is the gas density, 
v = (v x ,v y , v z ) T is the gas velocity vector, 

e = e+hv\ 2 = h e - V - (19) 
2 p 

is the specific total energy, h e is the enthalpy, and the 
internal energy e is related to the gas pressure by 

p = ( 1 -l)pe, (20) 

where 7 (=5/3) is the adiabatic index. 

The equations are solved by a finite volume scheme with 
operator splitting, which is based on the AUSMDV de- 
scribed by Wada & Liou (1997). Liou & Steffen (1993) 
developed a remarkably simple upwind flux vector split- 
ting scheme called 'advection upstream splitting method' 
(AUSM). It treats the convective and pressure terms of 
the flux function separately. The AUSMDV has a blend- 
ing form of AUSM and flux difference, and improves the 
robustness of AUSM in dealing with the collision of strong 
shocks. We extended it to third-order spatial accuracy us- 
ing MUSCL (van Leer 1977) with a total variation dimin- 
ishing limiter proposed by Arora & Roe (1997). Since this 
scheme has a great advantage due to the reduction of nu- 
merical viscosity, fluid interfaces are sharply preserved and 
small-scale features can be resolved as in the 'piecewise 
parabolic method' (PPM) of Woodward & Colella (1984). 
The AUSMDV scheme is, however, simpler and has a lower 
computational costs than PPM. The code has passed suc- 
cessfully several tests and handles very well weak and high 
Mach number shocks. In the Appendix we present the re- 
sults for two such tests, a standard one-dimensional shock 
tube and a point explosion in three dimensions. 

One of the main purposes of our study is to determine 
how far from the production sites can the product of stellar 
nucleosynthesis propagate into the IGM. Thus, we have to 
deal with very different length scales in our simulation. If 
we assume, e.g., a minimum resolution length comparable 
to the scale of an individual SN remnant, ~ 30 pc, and set 
the size of the simulation box to forty times the virial ra- 
dius of the halo, i.e. 60 kpc for h = 0.5, then about 2,000 
zones are needed per dimension. Even with massive su- 
percomputers, it is difficult to carry out such a simulation 
in three dimensions. We have therefore adopted a 3-D 
'nested grid method'. Our scheme is similar to Tomisaka's 
(1996) two-dimensional version: the general algorithm is 
based on the works of Berger & Oliger (1984) and Berger 
& Colella (1989). Six levels of fixed Cartesian grids were 
used, with every fine grid being completely covered by a 
coarser one. We named the grids LI (the finest), L2,...., 
L6 (the coarsest), as shown in the right panel of Figure |[ 
The box size of grid LI was set equal to 2r v i r , and the mesh 
size of the Ln grid to twice that of the Ln — 1 level. All 
grids were centered within each other, with the finest cov- 
ering the whole galaxy halo. Since the cell number is the 
same (128 x 128 x 128) for every level Ln (n = 1,2, ...,6), 
the minimum resolved scale is about 22 pc and the size of 
the coarsest grid is 96 kpc. Thus, the scheme has a wide 
dynamic range in the space dimension. 



The grids are connected by the transfer of conserved 
variables: the mass density p, the components of the 
momentum density pv, and the total energy density pe. 
Since each coarse cell is resolved by exactly 2 3 fine cells, 
the part of the coarse grid covered by the finer grid is 
overwritten with the arithmetic averages of these vari- 
ables on 8 fine cells. The boundary conditions for a fine 
grid are determined instead by monotonic interpolation 
of physical variables on the coarser grid. For the coars- 
est level L6 we adopt outflowing boundary conditions by 
imposing for each variable a zero gradient. Since the 
Courant-Friedrichs-Lewy condition requires the time-step 
of a coarser grid to be twice as long as the time-step of the 
next finer level, the finer grids must be calculated much 
more often. The finest grid (LI) is calculated first using an 
appropriate time-step At. After two time-steps of finest 
grid, the grid one level coaser is advanced by the time-step 
of 2Ai. The scheme is repeated recursively until all grids 
are advanced. Every time the sequence changes from a 
fine to a coarse grid, the data on the fine grid are averaged 
onto the coarse one. 

The center of the fixed NFW dark matter potential is 
located at the center of the simulation boxes, and the halo 
initial gas distribution follows equation (0). The halo gas 
is assumed to be non self-gravitating, of primordial com- 
position, optically thin, and in collisional ionization equi- 
librium. We assume that the metal ejecta either do not 
mix with the ambient gas or do so at late times (cf. Paper 
I). Even if some mixing occurs inside the hot cavity gas 
thereby increasing its metallicity, because of the very low 
density of this gas the cooling rate should remain much 
lower than in the outer cooling shell. The external IGM 
is at T = 10, 000 K, as expected from photoheating by 
the SN progenitors. We use the radiative cooling func- 
tion of Sutherland & Dopita (1993) for primordial gas, 
and include the effect of Compton cooling off microwave 
background photons, 

A c = (5.4 x 10~ 36 ergcm~ 3 s _1 ) XP (1 + z ) 4 T. (21) 

Here \ is the ionized fraction, z (=9) is the redshift, and 
T is the gas temperature. The cooling and heating terms 
are integrated implicitly. 

Mechanical feedback from SNe is a critical process in 
galaxy formation studies and simulations, as it modifies 
the composition and thermodynamical state of the am- 
bient gas. Different approximate prescriptions have been 
adopted in the past, largely imposed by numerical limita- 
tions and a poor appreciation of how to implement feed- 
back in hydrodynamical simulations (Katz 1992; Navarro 
& White 1993; Mihos & Hernquist 1994; Mori et al. 1997; 
MacLow & Ferrara 1999; Mori, Yoshii, & Nomoto 1999). 
Navarro & White (1993) and Mihos & Hernquist (1994) 
assumed that a fraction f v of the available SN explosion 
energy is deposited in the ambient g radial velocity 

perturbation directed away from the event, the remainder 
being dumped as heat. This method is obviously rather 
sensitive to the assumed value for /„. Mori et al. (1997) 
and Mori, Yoshii, & Nomoto (1999) assumed that the gas 
within the maximum radius of the shock front in the adi- 
abatic phase of a supernova remnant (Shull & Silk 1979) 
remains adiabatic until the multiple SNe II phase ends at 
t ms (8M Q ). In this case, the effects of radiative cooling 
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might have been be underestimated. MacLow & Ferrara 
(1999) modeled SN feedback in a central region of a dwarf 
galaxy as a constant luminosity wind driven by a ther- 
mal energy source. Nearby starbursts, however, are known 
to have multiple star formation sites scattered across the 
whole galaxy. 

Our algorithm for simulating SN feedback improves 
upon previous treatments in several ways. OB associa- 
tions are distributed as a function of gas density (oc p a ) 
using a Monte-Carlo procedure (the projected distribu- 
tions of associations in plotted in Figure . After a main 
sequence lifetime, all stars more massive than 8 M Q ex- 
plode instantaneously injecting an energy of 10 51 erg, and 
their outer layers are blown out leaving a compact rem- 
nant of 1.4 M Q . Therefore SNe inject energy (assumed in 
pure thermal form) and mass into the interstellar medium: 
these are supplied to a sphere of radius 

i^pcf^fV^^r 175 , (22) 
v V 10 y r / VO.lcm- 3 / V ; 

corresponding to the radius of a SN remnant in a uni- 
form ambient medium of density n and in the adiabatic 
Sedov- Taylor phase; the expansion velocity R s is also self- 
consistently calculated from such solution. The gas then 
starts cooling immediately according to the adopted gas 
cooling function. The radius R s has a minimum value of 
22 pc due to numerical resolution. The time-step At is 
controlled by the Courant condition. One additional con- 
straint is that, if there are more than two SN events in 
a OB association, we decrease At until only one explo- 
sion per association occurs during the time-step At. The 
main limitation of our method is the assumption that the 
gas is initially in hydrostatic equilibrium and non self- 
gravitating. 

5. NUMERICAL RESULTS 

As mentioned above, we have run two numerical sim- 
ulations, an extended stellar distribution (Case 1) and a 
more concentrated one (Case 2). In both runs 11170 SNe 
explode in our halo (corresponding to /* = 10%). Snap- 
shots of the gas density and temperature distributions in 
a slice along the X-Y plane of the nested grids are shown 
in Figures [| and ^| for Case 1, and Figures [?] and ^ for 
Case 2. The three columns in each figure depict the time 
evolution from about 5 Myr to up to 200 Myr. Along a 
given row, the leftmost panel refers to grid L5 (linear size 
48 kpc), the central one to grid L3 (linear size 12 kpc), 
and the rightmost panel refers to the grid LI (linear size 
3 kpc). The density range is — 5 < log (n/cm~ 3 ) < 1, and 
the temperature range is 4 < log (T/K) < 8.5. 

After a few Myr from the beginning of the simulation the 
most massive stars explode as SNe and produce expanding 
hot bubbles surrounded by a cooling dense (n « 1 cm~ 3 ) 
shell. At these early times the typical bubble size tends 
to be larger in Case 1 than in Case 2 because in the 
former scenario the explosions sites are more uniformly 
distributed and more likely to occur in the outer, lower- 
density regions, which favors their rapid expansion. In 
spite of the smaller typical bubble size, the degree of over- 
lapping appears to be more pronounced in Case 2, indi- 
cating that the crowding effect is dominant. 



The different initial topology of the multiphase ISM 
leaves an inprint also in the later evolutionary stages. At 
around 10 Myr the ISM structure is relatively ordered in 
Case 2, where individual bubbles have merged in a coher- 
ent (although far from spherical) expanding superbubble 
from which cold (T s» 10 4 K) filaments protrude inside 
the cavity. These filaments are the leftovers of previous 
individual shell-shell interactions further processed by hy- 
drodynamic instabilities. In Case 1, the halo topology is 
more perturbed, with bubbles expanding in the outer re- 
gions having already undergone blowout and venting their 
hot gas into the IGM (see, e.g., the structure at the top 
of the panel in the second row, third column of Fig. [|), 
and others whose interaction is giving rise to an intricated, 
multiphase structure in the inner halo, where 10 8 K gas co- 
exists with a cooler 10 4 K phase from which it is separated 
by cooling interfaces. Note that for Case 1 SN explosions 
in the outer halo drive inward-propagating shocks that act 
to collect and pile up gas towards the center. This effect is 
much less pronounced in Case 2 where the net mass flow 
is an expanding one. As we will see below, the impact of 
the mechanical energy deposition on the host pregalactic 
halo is rather different in the two simulations. 

As the evolution continues, a coherent and increas- 
ingly spherical shell expanding into the IGM is eventually 
formed in both runs. The shell contains a large fraction of 
the halo gas that has been swept-up during the evolution; 
at t = 20 Myr its size is roughly 6 kpc for Case 1 and 
slightly larger for Case 2, as a consequence of the most 
efficient use of mechanical energy in the latter simulation. 
In Case 1 one can clearly see a central, dense core resulting 
from the 'implosion' wave mentioned above. While such a 
feature is almost absent in Case 2 at t — 35 Myr, a dense 
core will form at later stages as a result of the accretion 
of cold clumps that are balistically raining towards the 
center. 

The final two bottom rows of the simulation figures show 
the final stages of the evolution that are qualitatively sim- 
ilar to the one just described. The shell is now nearly 
spherically symmetric: its interior is filled with warm 
(T 10 6 K) gas at a very low density n ^ 10~ 4 cm" 3 , 
i.e. slightly below the mean value for the IGM. Figure KJ 
shows the locus of the spherically-averaged shell radius 
r s heii as a function of time. The shell initially follows 
an energy-driven phase where r s h c ii oc t 3 / 5 (Weaver et 
al. 1977). At later time the evolution relaxes to the adia- 
batic Sedov-Taylor solution with r s h c ii oc t 2 ^ 5 . Afterwards 
the shell slows down to a momentum-conserving 'snow- 
plough' phase, r s heU oc i 1 / 4 . Figure || shows that, at the 
end of the simulations, the shell is still sweeping out IGM 
material; its radius and velocity are 21 kpc and 26 km 
sec -1 at t = 250 Myr. At Mach number M. = 1, as the 
shock will decay into sound waves, the maximum 'stalling 
radius' will be reached. Using momentum conservation we 
estimate the final radius of the shell to be close to 26 kpc. 
This is about a factor 2.3 larger than the value we derived 
analytically in Paper I for the same case. 

In both simulations runs the final configuration includes 
a central core resulting from the two different mechanisms 
already outlined; its central density and radial profile is 
similar to the initial one. However, the relative ratio be- 
tween the gas mass contained in the shell and the one in 
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the central core is different. A visual inspection of Figures 
H and ^ already shows a thicker shell and a less massive 
core in Case 2 relative to Case 1. This is confirmed quan- 
titatively by the plots in Figure [l(] (Case 1) and Figure 
|l2| (Case 2) where we show the fraction of the initial halo 
baryonic mass contained inside (1, 0.5, 0.1)r V i r as function 
of time. The differences are striking: in Case 1, the amount 
of gas at the center is constantly increasing, finally collect- 
ing inside 0.1 r V j r about 30% of the total initial mass. On 
the contrary, in Case 2 the central regions remain prac- 
tically devoided of gas until 60 Myr, when the accretion 
process starts. The final result is a small core containing 
a mass fraction of only 5%. In the former case 50% of the 
halo gas mass is ejected together with the shell, whereas in 
Case 2 this fraction is ~ 85%, i.e. the blow-away is nearly 
complete. 

The gas thermal history is illustrated in Figures [ll] (Case 
1) and [13| (Case 2), where the evolution of the filling factor 
of the components with temperature larger than a given 
threshold [log(T/K) = 4,5,6,7,8] is shown inside r v ; r . In 
both simulations we observe an initial conversion of cold 
gas into hot gas followed by the opposite, much slower, 
process. The main difference between the two runs con- 
sists in the larger extension of the warm T = 10 5 " 6 K 
component filling 20-50% of the volume in Case 2 at 120 
Myr. Also the compression suffered in the inner regions 
results for Case 1 in a very compact cold gas core, con- 
fined to 25% of the volume by 60 Myr, which re-expands 
afterwards. Case 2 produces a more extended hot gas dis- 
tribution as seen from the final values of the curves at 120 
Myr. 

Figure |lj shows the kinetic energy flux, , carried by 
the outflowing gas through the virial radius. The kinetic 
energy flux is calculated as 



4nr ■ F 



with 



F = J2 dV ilpiVr 3 i /J2 dVi > 



(23) 



(24) 



where dVt is a volume of the grid cell located at r v - lr from 
the center, and pi and v ri are the gas density and the radial 
velocity at that location, respectively. The summation is 
taken only for the case of a positive radial velocity at the 
virial radius. The mean SN mechanical luminosity, Lsn, 
is defined as 



Esn 



i ms (8Af ) - t u 



1.1 x 10 40 ergs 



5 (120M Q ) 
-l 



(25) 
(26) 



From Figure 14, we see that kinetic energy is ejected from 
the galaxy at roughly a constant fraction 25% of the SN 
mechanical luminosity from 13 Myr to 35 Myr for the Case 
1 run; this value increases to about 40% (but with a larger 
time spread) from 13 Myr to 35 Myr for the Case 2 run. 
Averaged over the entire evolution, we find that 23% (30%) 
of the total SNe energy is converted to the kinetic energy 
of the outflowing gas for the run of Case 1 (Case 2). 

6 We note that, on this assumption, and hence the early luminosity of these systems, would exceed the value predicted by the usual 
low-mass extrapolation of CDM galaxy formation models (cf. White & Frenk 1991). In these scenarios /* is postulated to decline steeply in 
shallow potential wells, thereby reducing the population of low— luminosity galaxies and avoiding the so— called 'cooling catastrophe'. 



6. SUMMARY AND DISCUSSION 

In this Paper we have used three-dimensional hydro- 
dynamic simulations to investigate the dynamics of SN- 
driven bubbles as they propagate through and escape the 
grasp of subgalactic halos with masses M = 10 8 h^ 1 M 
at rcdshift z — 9. Depending on the stellar distribution, 
we found that less than 30% of the available SN energy 
is converted into kinetic energy of the blown away ma- 
terial, and that mechanical feedback is less efficient than 
expected from simple energetic arguments, as off-nuclear 
SN explosions drive inward-propagating shocks that tend 
to collect and pile up cold gas in the central regions of the 
host halo. For the more extended star formation Case 1, 
the implosion wave collects more material in the center, 
eventually resulting in a more massive cold core in which 
star formation can be restarted. The amount of kinetic 
energy ejected into the IGM is essentially double for Case 
2, again because less energy flows towards the center in 
this run. Blow-away is more efficient if stars form prefer- 
entially at the center of the halo. 

SN-driven pregalactic outflows may be an efficient 
mechanism for distributing the product of stellar nucle- 
osynthesis over large cosmological volumes prior to the 
reionization epoch. In Paper I we used several approx- 
imations to show that, for a star formation efficiency of 
/* = 10%, the radius of the SN-driven bubble around our 
prototypical early halo would grow up with time up to a 
final stalling value of about 11 kpc, when pressure equi- 
librium with the IGM was achieved. In this simulations 
we find larger final radii, up to about 26 kpc for both 
runs. What is the expected spatial extent of the ensem- 
ble of wind-driven ejecta from a population of pregalactic 
systems? Consider an Einstein-de Sitter cosmology with 
Q M = 1, h = 0.5, cr 8 = 0.63, n= 1, and n b h 2 = 0.019, 
where the amplitude of the power spectrum has been 
fixed in order to reproduce the observed abundance of 
rich galaxy clusters in the local universe (we have used 
the transfer function by Efstathiou, Bond, & White 1992 
with shape parameter T = 0.5). According to the Press- 
Schechter formalism, the comoving abundance of collapsed 
dark halos with masses M w 10 s h~ 1 M ( r ) at z — 9 is 
then close to 80 ft. 3 Mpc -3 , corresponding to a mean proper 
distance between neighboring halos of 15 ft -1 kpc, a total 
mass density parameter of order 3 percent, and a stellar 
density parameter of 0.002 A ACDM cosmology with 
the same normalization would give a halo number density 
a factor 3.7 lower, resulting in a larger mean separation 
of 23 hr 1 kpc. Therefore, while with /* = 10% only a 
small fraction, about 4% percent, of the total stellar mass 
inferred at the present epoch (Fukugita, Hogan, & Peebles 
1998) would actually form at these early epochs, the im- 
pact of such an era of pregalactic outflows could be quite 
significant, as the product of stellar nucleosynthesis would 
be distributed over distance that are comparable to the 
mean proper distance between neighboring low-mass sys- 
tems, i.e. volume filling factors of order unity. 6 The col- 
lective explosive output of about ten thousands SNe per 
M ^ 10 s ft -1 M Q halo at these early epochs could then 
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pollute the entire intergalactic space to a mean metallicity 
(Z) = Qz/^b ft 0.003 (comparable to the levels observed 
in the Lya forest at z « 3) without much perturbing the 
IGM hydrodynamically, i.e. producing large variations of 
the baryons relative to the dark matter. 

We want to comment here on the possible effect of gas 
self-gravity, neglected in our simulations. The local gas- 
dynamical time of cold filaments and blobs in the halo 
ISM is of order 3 x 10 7 yr for a gas number density of 1 
cm~ 3 , and 9 x 10 6 yr for a density ten times greater. This 
timescale is comparable to the characteristic timescale of 
gas removal from the potential well of the dark matter halo 
in our simulations. If the density of the cold component is 
increased due to self-gravity, radiative cooling will be en- 
hanced. This will decrease the efficiency of the conversion 
of the available SN energy into the kinetic energy of the 
outflowing material, weakening the blow-away. Two com- 
pensating effects - that act to prevent gas cooling - should 
also be considered then, photoionization by UV radiation 
and thermal conduction between the cold and hot inter- 
stellar medium. In particular, the conduction timescale 
is 

3nkT 

Tcond = 2V • (kT 5 / 2 VT) 

- 4 - 3Myr (l^) (Kxb) 2 (nfK)" /2 ' (27) 

where k is the conduction coefficient, n is the gas num- 
ber density, and I is the characteristic scale length of 



the cold filaments. Taking T = 10 7 K, I = 100 pc, and 
n = 1 cm from the results of our simulations, we find 
Tcond = 4.3 x 10 6 yr. This timescale is shorter than the 
free-fall time and comparable to the cooling timescale. 
The thermal energy of the cold phase will then be affected 
by thermal conduction from the hot medium. 

Finally, while in this paper we simulated the effect of an 
instantaneous burst of star formation in a protogalactic 
spherical halo embedded in a uniform IGM, in a series of 
forthcoming studies we will investigate non-instantaneous 
stellar feedback in continuous star forming halos of differ- 
ent masses and morphologies, within a realistic cosmolog- 
ical density field. In these cases the blow-away history of 
the gas may be quite different. 
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APPENDIX 

A number of different tests were performed to demonstrate the ability of our code to reproduce known analytical 
results. Here we present results for: (i) a standard one-dimensional shock tube, and (ii) the adiabatic expansion of a 
point explosion in three dimensions. The hydrodynamic equations ( |l3| ) were solved with the source term S set equal to 
zero, on a grid of 128 zones in each dimension. 

In the shock tube two regions of different gas densities are instantaneously brought into contact (see Sod 1978). A 
shock (rarefaction) wave propagates then into the low (high) density gas. We choose the initial conditions to be: 

p=l, P=l, v = 0, for x < 0; , 

(0 = 0.125, P = 0.1, v = 0, for x > 0. [ ' 

The ratio of specific heats is 7 = 1.4 and the shock Mach number is 1.66. In this test we integrate the hydrodynamic 
equations on a fix ed g rid, i.e. without the nested grid method. The numerical results for the density, pressure, and velocity 
are shown in Fig. [L5| and compared with the analytical solutions. The code can resolve a shock front within 2 cells, while 
contact discontinuities are spread over 3 cells. Note how the AUSMDV scheme suppresses oscillations in the postshock 
velocity profile (cf. Cen 1992). 

The second test run follows the evolution of a spherical blast wave caused by a point explosion in a homogeneous gas 
with 7 = 5/3, neglecting radiative cooling. Two grid levels (LI and L2) are used in this case. The gas density is set 
to pi = 10~ 24 gr cm~ 3 and the gas temperature to 7$ = 100 K. An explosion energy of Ei = 10 51 ergs is added to the 
eight zones at the center of the LI grid. The propagation of the shock front is described by the self-similar Sedov-Taylor 
solution r s hcii = £o(i?»/ Pi)t 2 ^ 5 , where the coefficient £0 is 1.15 for 7 = 5/3. Fig. |l6| shows the normalized density profile as 
a function of radius at elapsed times 2638 yr and 10382 yr. The shock wave propagates outward from the center of the LI 
grid, and smoothly passes through the boundary between grid levels LI and L2. The numerical results nicely reproduce 
the analytic shock position and profile. 




Fig. 1. — Isocontours in the xq TV piano of bubble shell velocity, v B (r v i T ) (in km s — 1 ), at the virial radius. Here N SNe are assumed to 

explode at a distance xo from the center (in units of the halo virial radius). Also shown (dotted line) is the isocontour corresponding to the 
escape speed v e (r v i T ). The points indicate the locus of the associations used in the numerical simulations for a = 1 (Case 1). 
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Fig. 3. — Halo gas column density S(xo) (normalized to the total gas column density through the halo St) vs. [fs(r v i r )/t!e(r v ir)] 2 for two 
values of TV = 10, 100 (solid curves). The dotted lines depict the corresponding curves of constant kinetic energy flux, £ k (N). The shaded 
area indicates the parameter space in which the wind is energy-carrying. 
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Fig. 4. — Le/i: Projected distributions of OB associations. The upper and lower panels correspond to Case 1 and Case 2, respectively. 
Right: The structure of the nested grids is shown for a two-dimensional sectional plane. Six levels of fixed Cartesin grids were used, with the 
same number (128 X 128 X 128) of cells for every level. The box size of the finest grid (LI) was set equal to 2r v i r , so LI covers the entire 
subgalactic halo. 
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Fig. 6. — Snapshots of the logarithmic temperature of the gas at five different elapsed times for our Case 1 simulation run. The three panels 
in each row show the spatial temperature distribution in the X — Y plane on the nested grids. The left, middle and right panels in each row 
correspond to the level L5, L3, and LI grid, respectively. 
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Fig. 7. — Same as in Fig. B but for our Case 2 simulation run. 
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9. — The evolution of the shell radius as a function of time. The filled (open) circles correspond to the Case 1 (Case 2) simulation 



run. The solid lines show the expected power-law behavior in the energy-driven phase (r she n oc t 3 / 5 ), the adiabatic Sedov-Taylor solution 
( r shell <x * 2//5 )> anc l the momentum-conserving 'snowplough' phase (r s h e ii oc t 1 / 4 ). 
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Fig. 10. — The evolution of the gas mass inside the gravitational potential well of the CDM halo for our simulation run of Case 1 (a = 1) as 
a function of time. Curves correspond to the gas mass inside the virial radius r v ; r (solid line), 0.5r v ; r (dashed line), and 0.1r v ; r (dotted line). 
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Fig. 11. — Cumulative evolution of the volume filling factor inside the virial radius of the gravitational potential well of the CDM halo for 
our simulation run of Case 1 (a = 1), as a function of time. Each curve corresponds to a different gas temperature (solid line: T = 10 4 K, 
dashed line: T = 10 5 K, dash-dotted line: T = 10 6 K, dotted line: T = 10 7 K, and dash-dot-dotted line: T = 10 8 K). 
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Fig. 14. — Fraction of the SN mechanical luminosity carried by the outflowing gas as kinetic energy through the virial radius, as a function 
of time. The thick line corresponds to the simulation run of Case 1, and thin line corresponds to that of Case 2. 
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Fig. 16. — Point explosion in a homogeneous medium. The normalized density profile is shown as a function of the (self-similar) radius at 
elapsed times of 2638 yr (filled circles) and 10382 yr (open circles). The self-similar Sedov-Taylor solution for the gas density is indicated by 
the solid line. The inlet shows a two-dimensional cut through the explosion site. 



